import std.random;


struct Pos
	{
	int x, y;
	}

struct Z
	{
	uint z1, z2, z3, z4;
	};

uint lcgStep(ref uint z, const uint A, const uint C)   
	{return z = A * z + C;}  
// S1, S2, S3, and M are all constants, and z is part of the  
// private per-thread generator state.  
uint tausStep(ref uint z, const int S1, const int S2, const int S3, const uint M)  
	{  
//	const uint b=(((z << S1) ^ z) >> S2);  
	return z = (((z & M) << S3) ^ (((z << S1) ^ z) >> S2));  
	}

uint hybridTausU(ref Z z)  
	{  
	// combined period is lcm(p1,p2,p3,p4)~ 2^121  
	return   
		tausStep(z.z1, 13, 19, 12, 4294967294UL) ^  // p1=2^31-1  
		tausStep(z.z2, 2, 25, 4, 4294967288UL) ^    // p2=2^30-1  
		tausStep(z.z3, 3, 11, 17, 4294967280UL) ^   // p3=2^28-1  
		lcgStep(z.z4, 1664525, 1013904223UL)        // p4=2^32  
		;  
	}  

const int nPop = 256;
const int nInd = 256;
const int maxX = 1024;
const int maxY = 1024;
const int cellSz = 64;
const int cacheMX = maxX / cellSz;
const int cacheMY = maxY / cellSz;
const int vis = 10;
const int maxSqDist = vis * vis;

int cachePos(int x)
	{
	return x / cellSz;
	}

int m2i(in int x, in int y)
	{
	return y * cacheMX + x;
	}


struct Cache(T)
	{
	T[] data;
	int n;

	void add(T p)
		{
		if (n >= data.length)
			data.length = data.length + 1;
		data[n] = p;
		n++;
		}

	void remove(T p)
		{
		for (int i=0; i<n; i++)
			{
			if (data[i] == p)
				{
				data[i] = data[n-1];
				n--;
				return;
				}
			}
		}
	}

int main()
	{
	Pos[][] pos;
	Cache!(int)[][] cache;
	int[][] coll;

	pos.length = nPop;
	cache.length = nPop;
	coll.length = nPop;

	for (int p=0; p<nPop; p++)
		{
		pos[p].length = nInd;
		coll[p].length = nInd;
		cache[p].length = cacheMX * cacheMY;
		for (int i=0; i<nInd; i++)
			{
			pos[p][i].x = uniform(0, maxX);
			pos[p][i].y = uniform(0, maxY);
			cache[p][m2i(cachePos(pos[p][i].x), cachePos(pos[p][i].y))].add(i);
			}
		}

	Z z;
	z.z1 = uniform(0, 30000);
	z.z2 = uniform(0, 30000);
	z.z3 = uniform(0, 30000);
	z.z4 = uniform(0, 30000);

	for (int s=0; s<1000; s++)
		{
		for (int p=0; p<nPop; p++)
			{
			for (int i=0; i<nInd; i++)
				{
				const int cpX = cachePos(pos[p][i].x);
				const int cpY = cachePos(pos[p][i].y);

				const int newX = (pos[p][i].x + hybridTausU(z)%5 - 2 + maxX) % maxX;
				const int newY = (pos[p][i].y + hybridTausU(z)%5 - 2 + maxY) % maxY;

				const int cpnX = cachePos(newX);
				const int cpnY = cachePos(newY);

				if (cpnX != cpX || cpnY != cpY)
					{
					cache[p][m2i(cpX, cpY)].remove(i);
					cache[p][m2i(cpnX, cpnY)].add(i);
					}
				
				pos[p][i].x = newX;
				pos[p][i].y = newY;
				}

			for (int i=0; i<nInd; i++)
				{
				const int minCpX = cachePos(pos[p][i].x - vis);
				const int minCpY = cachePos(pos[p][i].y - vis);
				const int maxCpX = cachePos(pos[p][i].x + vis);
				const int maxCpY = cachePos(pos[p][i].y + vis);

				for (int cpx=minCpX; cpx<=maxCpX; cpx++)
					{
					const int offs_x = cpx<0 ? -maxX : 
						(cpx>=cacheMX ? maxX : 0);
					const int rx = (cpx+cacheMX)%cacheMX;
					for (int cpy=minCpY; cpy<=maxCpY; cpy++)
						{
						const int offs_y = cpy<0 ? -maxY : 
							(cpy>=cacheMY ? maxY : 0);
						const Cache!(int) * local = &(cache[p][m2i(rx, (cpy+cacheMY)%cacheMY)]);
						for (int l=0; l<local.n; l++)
							{
							const int dx = pos[p][i].x - (pos[p][local.data[l]].x + offs_x);
							const int dy = pos[p][i].y - (pos[p][local.data[l]].y + offs_y);
							const int sqDist = dx*dx + dy*dy;
							if (sqDist < maxSqDist)
								coll[p][i]++;
							}
						}
					}
				}
			}
		}


	return 0;
	}
